Lattice dynamics of mixed semiconductors (Be,Zn)Se from first-principles calculations 
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Vibration properties of Zni-^^BeajSe, a mixed II- VI semiconductor characterized by a high con- 
trast in elastic properties of its pure constituents, ZnSe and BeSe, are simulated by first-principles 
calculations of electronic structure, lattice relaxation and frozen phonons. The calculations within 
the local density approximation has been done with the Siesta method, using norm-conserving 
pseudopotentials and localized basis functions; the benchmark calculations for pure endsystems 
were moreover done also by all-electron WIEN2k code. An immediate motivation for the study was 
to analyze, at the microscopic level, the appearance of anomalous phonon modes early detected 
in Raman spectra in the intermediate region (20 to 80%) of ZnBe concentration. This was early 
discussed on the basis of a percolation phenomenon, i.e., the result of the formation of wall-to-wall 
-Be-Se- chains throughout the crystal. The presence of such chains was explicitly allowed in our 
simulation and indeed brought about a softening and splitting off of particular modes, in accor- 
dance with experimental observation, due to a relative elongation of Be-Se bonds along the chain 
as compared to those involving isolated Be atoms. The variation of force constants with interatomic 
distances shows common trends in relative independence on the short-range order. 

PACS numbers: 63.20.-e, 63.50.-l-x, 71.15.-m, 71.55.Gs, 78.30.-j 
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I. INTRODUCTION 



(Be,Zn)Se is an example of mixed II-VI semicon- 
ductor system whose electronic properties and, conse- 
quently, elastic characteristics vary with concentration 
in a non-trivial way. A considerable effort has been 
spent on characterizing the optical band gap, which drops 
slightly from intermediate concentrations and undergoes 
a change from direct to indirect character underway from 
ZnSe to BeSei. Not less interesting possibilities open in 
the tuning of elastic properties, because BeSe and ZnSe, 
similarly to some other easily miscible II-VI constituents, 
possess quite different elastic constants, which have to 
be accommodated in a solid solution. (Be,Zn)Se makes, 
actually, quite an extreme case of elastic contrast be- 
tween constituents in a mixed system. This was argued 
to favor the appearance of anomalous lines in the Ra- 
man spectra, which split off on the soft side of the BeSe- 
related TO peak in the x ~0.2-0.7 concentration range 
of Bea;Zni_j;Se crystals, as was detected by Pages at al? . 
In order to explain this anomaly, Pages et al? brought 
into discussion an idea of quasi-infinite chains, which are 
developing in the softer matrix (ZnSe) as the BeSe con- 
centration reaches the percolation threshold of ~20%. 
However, the microscopic mechanism relating percola- 
tion with vibration properties remained speculative. The 
present work's aim is to provide a possible complete and 
reliable first-principles description of internal tensions, 
structural relaxation, and the impact of the latter on 
lattice-dynamical properties of a mixed Bea;Zni_a;Se al- 
loy near the percolation threshold. A concise result of 
our simulation has been reported earlier as a conference 



proceeding^. Here we offer a detailed outline of results 
obtained for different supercells, with a more lengthy dis- 
cussion. In particular, we provide an analyzis of bond 
lengths distribution for a number of compositions, dis- 
cuss the change of force constants related with that, and 
finally calulate phonon density of states and provide its 
the wavevector-resolved decomposition. 

The electronic properties of pure constituent com- 
pounds are well known; they also were probed by ab 
initio methods of respectful accuracy in a number of ear- 
lier publications. Simulations of comparable accuracy on 
mixed alloys are seriously complicated by the necessity 
to treat large supercell, with a practical loss of any useful 
symmetry. A certain success has been achieved in molec- 
ular dynamics (MD) simulations of an isostructural III-V 
mixed semiconductor alloy, (Ga,In)As: Branicio et al^ 
performed parametrized MD calculations for supercells 
with up to 8000 atoms, and ab initio MD calculations on 
216-atom supercells. We are not aware of any elastic sim- 
ulations on mixed II-VI semiconductors but that by Tsai 
et al&, who used multicenter MD method, a formally ab 
initio one but of inferior accuracy to the method applied 
here. Comparing our results with those of Tsai et al., we 
find an overall good agreement, albeit with quantitative 
discrepancies on some sensitive issues. A more important 
difference is, however, that we study more attentively the 
limit of low Be concentration, particularly the onset of 
percolation on the Be sublattice. Moreover, we calculate 
phonon densities of states and discuss them in immediate 
reference to experimentally obtained Raman spectra. 

The paper is organized as follows. In Sec.|n]we outline 
our ab initio calculation approach to electronic structure 
and vibration properties. Sec. IIIII containing the results 
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for pure endsystems ZnSe and BeSe, is concise in view 
of numerous previous studies, but it is necessary here in 
order to assess the accuracy of the calculation method. 
In Sec. IIVI we outline our results for lattice parameters 
and equilibrium bond lengths in relaxed supercells, sim- 
ulating a broad concentration range. The effect of bond 
lengths on the force constants between nearest and next- 
nearest neighbors is discussed in Sec. Finally, the 
manifestation of the force constants in phonon frequen- 
cies is explained in Sec. IVll along with the discussion on 
wavevector dependency and vibration patterns in differ- 
ent parts of the phonon spectrum. 



II. CALCULATION SCHEME AND SETUP 

Dynamical properties are derived from the total en- 
ergy or forces, which are evaluated ab initio in a 
sequence of density functional calculations. We ap- 
plied the calculation method, and computer code, 
SiESTi^, which incorporates norm-conserving pseudopo- 
tentials in combination with atom-centered strictly 
confined numerical basis functionsSiS. The pseu- 
dopotentials were constructed along the TrouUier- 
Martins schemeiS for the following valence-states con- 
figurations: Be2s2(2.0), Se3di°(1.2) 432(1.9)4^4(2.0), 
Zn3di"(1.09)4s2(2.28). The numbers in brackets stand 
for pseudoization radii of corresponding states in Bohr. 
The basis set consisted of double-^ functions with polar- 
ization orbitals in all I channels on all centers (Be2s2p, 
Zn3(i4s4p, Se3d4s4p), hence also for exphcitly included 
semicore Se3d states. The details on the basis classifica- 
tion and testing in Siesta can be found in Refs.|M9. Our 
basis generation followed the standard split schema of the 
Siesta code, with the energy shift parameter, responsible 
for the localization of basis functions, equal to 20 mRy. 
This resulted in basis functions with maximal extension 
of 3.17 A (Be), 3.12 A (Zn) and 2.82 A (Se), i.e. weU be- 
yond the nearest-neighbor distances but short of direct 
overlap to next-nearest neighbors. The full structure op- 
timization (of cell parameters and internal coordinates) 
was performed in all supercells prior to lattice dynamics 
calculations. In order to get reliable forces on atoms it is 
essential to accurately perform spatial integration of the 
(general-form) residual charge density over the unit cell. 
To this end, we used the mesh cutoff parameter of 250 Ry, 
that generated a real-space mesh with the step of about 
0.1 A along each Cartesian direction throughout the su- 
percell. Within each mini-cell of this real-space mesh, 
an averaging of integration results was done over four 
fcc-type sampling points for better numerical stability. 
This resulted in the forces summing up, over all atoms in 
the supercell, to zero within the accuracy of ±0.1 eV/A, 
separately along each Cartesian direction. The thus vali- 
dated forces were used to calculate the dynamical matrix 
elements, by introducing small (0.016 A) deviations of 
atoms one by one from their equilibrium positions. The 
details of phonon calculations are discussed in Sec. IVII 



All our Siesta calculations have been done in the lo- 
cal density approximation (LDA). For the pure ZnSe and 
BeSe systems, moreover, we used the full-potential aug- 
mented plane wave method (see, e.g., Ref. ITlil . as im- 
plemented in the WIEN2k code^^. We apply here this 
all-electron method of recognized accuracy for bench- 
mark calculations of elastic properties. We provide the 
WIEN2k results obtained both within the LDA and in 
the generalized gradient approximation (GGA, Ref. 113) . 



III. ENDSYSTEMS, ZnSe AND BeSe 

The electronic structure of zincblende-structure semi- 
conductors ZnSe and BeSe was extensively studied by 
both experimental means and first-principle calculations. 
Band dispersions and partial densities of states have been 
calculated in the LDA by Lee et alm^ and Agrawal et 
alM> (for ZnSe, among other Zn chalcogenides) and by 
Fleszar and Hankei^ (for BeSe, among other Be chalco- 
genides). In the latter publication, also the quasiparticle 
band structure has been calculated in the GW approx- 
imation, and subsequently discussed in the context of 
LDA vs. GW and exact exchange for a number of sp 
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FIG. 1; The total energy vs. lattice constant curves (left) 
and the energy profiles corresponding to the F-TO phonon 
(right) as calculated for ZnSe and BeSe. A^, indicates the 
magnitude of the cation (or anion) displacement along [100] 
from the equilibrium. Open squares: WIEN2k (LDA) results; 
connected black dots: Siesta (LDA) results. 
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TABLE I: Elastic properties of ZnSe and BeSe from experiments and first-principles calculations. 







ZnSe 






BeSe 




Method 


a (A) 


B (Kbar) 


uto (cm ^) 


a (A) 


B (Kbar) 


ujto (cm~^) 


WIEN2k (LDA) 


5.568 


714 


198 


5.084 


920 


498 


WIEN2k (GGA) 


5.571 


727 


206 


5.182 


816 


523 


Siesta (LDA) 


5.587 


758 


203 


5.078 


912 


500 


exp. 


5.668" 


624'' 


207,' 205"^ 


5.137'= 


920' 


501'* 


other calc. 


5.677/ 5.638,9 
5.633/ 5.636*^ 


689/ 652,9 
811,'' 649*^ 


224'^ 


5.037' 


988' 


547^' 



"Ref. 
''Ref. 
'^Ref. 
'^Ref. 
^^Ref.] 

•''LDA, no rm- conserving pseudopotential with Zn3d as valence 
states, Ref.0 

'LDA, norm-cons ervi ng pseudopotential with Zn3d and Se 3d as 

valence states, Ref. |14| 

''LDA, full-potential LMTO, Refill 
'LDA, norm-conserving pseudopotential, 
■'LDA, norm-conserving pseudopotential, 
*^fuU-potential LMTO, Ref ^2^ 



zincblende seiniconductor32&. For ZnSe, the quasiparti- 
cle band structure has been reported by Luo et alm^. Our 
calculations provide the band structures in agreement 
with good earlier LDA results, so we skip the discussion 
on this point. Instead we turn to elastic properties and 
note that the optimized ground-state volume along with 
the bulk modulus have been reported in many calcula- 
tions for ZnSe (see Table 1), and also by Muhoz et alm^ 
and Gonzales-Diaz et alm^ for BeSe. The latter calcula- 
tions largely overestimate the stiffness of the BeSe crys- 
tal, probably due to the attribution of the Se3c? states to 
the core. As the treatment of Zn3d (and also preferably 
Se3c?) as valence states is now recognized as essential in 
pseudopotential calculations for achieving good descrip- 
tion of elastic properties of ZnSe, an overall agreement is 
established between state-of-art LDA calculations. Fig.Q] 
shows the total energy profiles, as function of uniform 
lattice scaling (left panels) and in its dependency on the 
off-center displacement (along any Cartesian direction) of 
one sublattice (anions or cations) relative to the other, 
i.e., the zone-center TO mode (right panels). The energy 
vs. lattice constant curves from all-electron WIEN2k cal- 
culations and from Siesta are almost identical. Small 
kinks in the Siesta total energy plot at some values 
of lattice constant are due to changes of the real-space 
mesh for spatial integration (as the cell volume varies 
but the mesh density is maintained about the same). 
The energy profiles of the anion-cation displacement (TO 
phonon, see right panels of Fig. ^) shows no substantial 
anharmonicity (deviations from parabolic behavior) for 
both ZnSe and BeSe. The phonon frequencies in Table |l] 
were calculated from the second-order fit to these data. 
Phonon frequencies in the following sections result from 
the diagonalization of dynamical matrices, constructed 



from the forces induced on atoms by small Cartesian dis- 
placements. For pure constituents they agree very well 
with the results of the second-order total energy fit. It is 
noteworthy that a big difference in F-TO frequencies of 
ZnSe and BeSe is only in part due to a large differences 
in the cation masses. The force constants are larger in 
BeSe, as is well seen from Fig. ^ indicating higher stiff- 
ness of this material, also manifested in its larger bulk 
modulus. 

Some our results listed in Table 1 are slightly refined 
as compared to those reported earlier in Ref. 4, due to an 
use of a more recent realization of the all-electron method 
(WIEN2k instead of WIEN97), and a more extended ba- 
sis set in Siesta. 

In the discussion of phonon properties below we re- 
fer to dispersion curves of ZnSe (obtained from neutron 
scattering and calculated in the rigid ion model, Ref. 
and BeSe (calculated using self-consistent pseudopoten- 
tial method and linear response approach, Ref. l29j) . A 
recent calculation^? of phonon dispersion for both sys- 
tems based on a simple model with central and angular 
forces yields acceptable results for ZnSe but large differ- 
ence from ab initio phonon dispersion curves for BeSe. 
This could be a manifestation of more covalent character 
of bonding in BeSe, as discussed below. 



IV. BOND LENGTHS IN MIXED CRYSTALS 

In a sequence of supercell calculations which represent 
mixed Be-Zn compositions, we begin by unconstrained 
relaxation of lattice vectors and internal coordinates. 
Probably the simplest model of a mixed system with pre- 
dominantly either Be or Zn at cation sites is an ordered 
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FIG. 2: Cubic supercells used in the calculation. The con- 
tinuous Be-Se chain and isolated Be atoms are marked in the 
large supercell. The figure was created with the XCrySDen 
software^--. 

8-atoni superstructure with the cubic primitive cell. Such 
supercells BeZn3Se4 and Be3ZnSe4 are shown in the left- 
hand side of Fig. ^ We emphasize that the BeZn3Se4 
superstructure contains Be atoms only in configuration 
which will be further referred to as "isolated" one; i.e., 
each Se neighbor to a Be atom has only Zn ions to sat- 
urate its other bonds. In the Be3ZnSe4 superstructure, 
on the contrary, each Se ion has but a single Zn neigh- 
bor, whereas Be-Se bonds form a continuous framework 
throughout the crystal. 

In order to better understand the changes of the 
phonon spectra in presumed dependency on the perco- 
lation setup, we performed the lattice dynamics simu- 
lation in a larger 2x2x2 (64 atoms) supercell. It con- 
tained 4 beryllium atoms in a continuous -Be-Se- chain 
transversing the crystal; moreover, two other Be atoms 
were situated in "isolated" positions, with only their 
fifth-nearest neighbors being of the Be type. This su- 
percell is also shown in Fig. |21 When fully relaxed, it 
maintains a nearly cubic shape, with a tiny z-compression 
(c/a=0.998) due to the anisotropy of chains. The nom- 
inal composition of Be is 18.75 at.%, i.e. just below 
the theoretical percolation threshold (0.198) on the fee 
lattice'^^; however, our structure includes percolation on 
the Be sublattice by construction. The advantage of our 
choice of supercell is that, despite its relatively moderate 
size, it allows to analyze local properties of areas with 
percolation (chained Be-Se bonds) or without percola- 
tion (isolated Be-Se bonds only) on equal footing. The 
results for this supercell will constitute the bulk of our 
discussion on the lattice dynamics. 

We note in passing that possible manifestation of per- 
colation effects in mixed semiconductor systems, notably 
at the above critical concentration and especially in sys- 
tems with large contrast in stiffness, has been pointed at 
by Bellaiche et al. already some time ago^. 

Finally, four supercells Be„Zn32-n Se32 (n=l,. . . ,4) of 
the same size and cubic shape as the latest one simulate 
a gradual aggregation of Be ions neighboring the same 



anion site. 

Fig. 01 shows the equilibrium lattice constant (derived 
from the relaxed supercell volume) and the mean cation- 
anion bond lengths. The mean lattice parameter exhibits 
a markedly linear dependence on the concentration. At 
the same time, the Zn-Se and Be-Se bond lengths tend 
to remain nearly constant throughout the whole concen- 
tration range. A similar conclusion follows from the MD 
simulations for the (Ga,In)As system by Branicio et aL-, 
who also cite EXAFS data, agreeing well with their re- 
sults, and for the (Be,Zn)Se alloy - by Tsai et al&, also 
based on MD simulations. However, the later publica- 
tion reports a seemingly too large drop in the Be-Se dis- 
tance when going from Zni/4Be3/4Se to pure BeSe. It 
should be noted that Tsai et al. found much smaller 
difference between the Zn-Se and Be-Se bond lengths 
than in our case, moreover their lattice parameters for 
pure constituents considerably deviate from experiment, 
probably, due to insufhciency of the basis set (and/or 
pseudopotential) they used. 

The diversity of local order in a mixed crystal com- 




BeSe Beo^^ZrigggSe BeoggZH^^gSe ZnSe 

FIG. 3; Top panel: calculated equilibrium lattice constant 
in Bea;Zni_a;Se crystals of different concentration, represented 
by supercells (see text for details). Bottom panel: equilibrium 
Be-Se and Zn-Se bond lengths. For x—0.19, the bond length 
values show a scattering due to the presence of different local 
environments - see details in the text and in Fig. U] 
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Be-Se bonds (A) 



2.4 

Zn-Se bonds (A) 



FIG. 4: Bond lengths in relaxed (Be,Zn)Se supercells; A 
discrete set of 128 nearest-neighbor distances is artificially 
smeared for better visibility. Top panels: bond lengths in the 
"model" supercell, with four Be atoms in continuous Be-Se 
chains, and two isolated Be impurities. Distances between 
difi'erent types of anions and cations are shown separately. 
Bottom panels: Be-Se and Zn-Se interatomic distances in 
four supercells, containing one to four Be atoms, neighboring 
the same Se atom. 



petes with the tendency to maintain the bond lengths, 
with the effect that the latter get a certain scattering 
around their mean values. This is illustrated in Fig. |2| 
for the Be6Zn26Se32 supercell, and in more detail for this 
and other "large" supercells - in Fig. 21 A discrete "spec- 
trum" of, in total, 128 bond lengths in each fully relaxed 
supercell is artificially broadened there for better visibil- 
ity. We start our discussion from the BeZn3iSe32 super- 
cell. Four Se neighbors experience an inward relaxation 
to a single Be impurity, shortening the Be-Se bonds to 
about 2.24 A, i.e., only slightly larger than in pure BeSe. 
Simultaneously the bonds to their Zn neighbors extend 
beyond the average interatomic distance (^2.41 A) in 
the rest of the predominantly-ZnSe supercell. By grad- 
ually adding more Be neighbors which share the same 
Se atom, we end up with a symmetric relaxation pat- 
tern in Be4Zn26Se32 with a set of 12 contracted bonds 
(between Be and their outer Se neighbors) and 4 ex- 



tended bonds (from each Be to the central Se). This 
is clearly seen in Fig. 0] The intermediate cases of 2 
and 3 Be atoms in the supercell offer an interpolation 
between the two discussed cases, introducing a diversity 
in the Be-Se bond lengths depending on the type of sec- 
ond neighbors. The splitting of the Be-Se bonds into 
shorter and longer ones, around the mean value of ~2.24 
A, leads to the diversification of the Zn-Se bonds into, 
correspondingly, longer and shorter ones, as compared to 
the dominating bulk-ZnSe central peak at ~2.41 A. A 
more complex distribution of bond lengths comes about 
in a supercell which contains, along with isolated Be im- 
purities, a continuous -Be-Se- chain (Fig.^ top panels). 
A Be atom substituting Zn favors an inward relaxation of 
neighboring Se. Large enough around a zero-dimensional 
defect (isolated atom), such inward relaxation is even 
stronger around a one-dimensional defect, such as an 
extended -Be-Se- chain. The tendency for shortening 
the Be-Se bond lengths should obviously exist also along 
the chain, but the actual contraction is limited by the 
fact that the chain is quasi-infinite and embedded into 
the ZnSe lattice, which effectively fixes the chain's step. 
However, since the repeated cation-anion sequence in the 
zincblende structure is folded, the shortening of its links 
can be achieved by accommodating the interbond angles, 
at the price of stretching the bonds between the in-chain 
Se atoms and their Zn neighbors in the crystal (shaded 
area in the top right panel of Fig. . The distribution of 
bond lengths between Zn cations and the off-chain Se an- 
ions resembles roughly that in the chain-free Be4Zn28Se32 
supercell. 



V. FORCE CONSTANTS 

We look now at how the variations in bond lengths 
map onto changes in the force constants. As already 
mentioned, the latter become accumulated as we displace 
the atoms one by one from their equilibrium positions 
along three Cartesian directions and analyze the forces 
induced on all atoms of the supercell. Specifically (in the 
symmetrized form). 
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a/3 



FtaR}+d^)-Fnm-d^) 



2d': 



Ff{{Il} + df)-F^^{{R}-df) 
2d? 



(1) 



where F"' is the force on atom a in the direction i, and 
{'R}+dj means that of all atoms, only the atom f3 is 
displaced along j from its equilibrium position, by (in our 
case) d=0.03 Bohr. With all elements of DfJ^ recovered, 
the solution of dynamical equation 
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FIG. 5: Diagonal elements of force constant matrix between 
nearest neighbors. See text for details. 



yields zone-center plionon frequencies lu and eigenvec- 
tors A" of the supercell. Before turning to the anal- 
ysis of phonons, we discuss briefly the force constants 
Dff ■ For each atomic pair {a, /3}, they make a 3x3 ma- 
trix which can be diagonalized. The diagonal elements 
are shown in Fig. [S] in dependency on the correspond- 
ing interatomic distance. In accord with previously dis- 
cussed scattering of the Be-Se and Zn-Se bond lengths, 
we have now a scattering of force constants. The major 
elements (of diagonalized 3x3 matrix) show a remark- 
ably pronounced linear variation with the bond length. 
Two minor elements are much smaller in magnitude and 
roughly bondlength independent. One can interpret the 
smallness of minor elements of D"^ as a measure of the 
covalency of corresponding bonds. Indeed, in case of a 
purely ionic bonding only the central motion of atoms 
would bring about a change in the force; the tangential 
component makes a non-zero contribution to the force 
constant only if tangential displacement induces a re- 
distribution of charge density, meaning that the cova- 
lency of the bond is not negligible. The contrast between 
strongly anisotropic and hence essentially ionic force con- 
stants for the Zn-Se bonds and much more covalent Be- 
Se interaction is clearly seen in Fig. [S] It is remark- 
able that the central Be-Se interaction is, on the average, 
smaller than the Zn-Se one. Yet, we have seen that the 
BeSe crystal possesses larger bulk modulus and a higher 
elasticity parameter in the F-TO vibration. This can 
only be explained by a considerable contribution of non- 
central forces (minor elements of the force constants) in 
the shaping of such elastic parameters. In other words, 
whereas in ZnSe the stretching of Zn-Se bonds plays a 
dominant role in the elasticity, in BeSe the bending of 
bond angles is nearly as important. 

We find very remarkable a pronounced linear depen- 
dence of the force constants on the interatomic distance, 
in apparent indifference to other aspects of the local or- 
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3.7 3.8 3.9 4 3.6 3.7 3.8 3.9 4 4.1 4.2 



Distance between next nearest neighbors (A) 

FIG. 6: Diagonal elements of force constant matrix between 
next nearest neighbors. Left panel: (Be,Zn)-(Be,Zn); right 
panel: Se-Se. 

der. The obtained dependence can be parametrized and 
used for fast lattice-dynamics simulations in much larger 
supercells. The results for similarly diagonalized second- 
nearest neighbors interaction are shown in Fig. The 
scattering of the Se-Se distances is twice as large as that 
of cation-cation ones, and shows a clear splitting into 
two groups, similar to that reported in Ref.y for the As- 
As pair distribution function in (Ga,In)As. This means 
that Be substitutes Zn almost at their original (non- 
displaced) positions, i.e., the cation sublattice remains 
rather rigid, whereas the Se atoms undergo large dis- 
placements, depending on their local environment. We 
see a pronouncedly different behavior of cation-cation 
and anion-anion force constants: the former have two 
almost symmetric (positive/negative) diagonal elements, 
with the third being nearly zero, independently on both 
type of pair (Be-Be, Be-Zn or Zn-Zn) and on the dis- 
tance. The anion-anion coupling has an appreciable 
central-force element, decreasing with the distance, and 
two much smaller tangential elements. The variation of 
the major element with distance is (at least in part) re- 
lated to the fact that on the left hand side of the plot 
we deal mostly with Se-Be-Se connections, in which the 
central Se-Se interaction reflects the effect of bending the 
strongly covalent Be-Se bonds; on the right hand site, the 
Se-Zn-Se connections offer a much weaker resistance to 
such bending, due to a higher bond ionicity. 



VI. PHONONS 

Since we deal with relatively large "disordered" su- 
percells with no internal (short-scale) periodicity, and 
no clear-cut phonon dispersion exists in such systems, 
it does not make sense to Fourier-transform the force 
constants ^ prior to solving the dynamical equation 
However, we explain below how to extract some 
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FIG. 7: Phonon density of states for the BeZn3Se4 and 
Be3ZnSe4 supercells, calculated for g=0 of the supercell and 
broadened by 10 cm~^. 



g-resolved information about lattice vibrations. The pri- 
mary characteristic for our discussion will be the phonon 
density of states (PhDOS), resolved when necessary over 
freely chosen groups H of atoms a: 



(3) 



This is a discrete spectrum of (37V — 3) lines (for TV atoms 
in the supercell, with acoustic modes removed) of differ- 
ent intensity, which is in the following figures broadened 
with the halfwidth parameter of 10 cm~^, for better vis- 
ibility. The vibration modes obtained from the solution 
of Eq. I^J correspond to the zone-center of the supercell 
in question, but they reflect different vibration patterns, 
also those of non-zone-center character, with respect to 
the underlying zincblendc lattice. Let us discuss this for 
8-atom supercells, whose calculated PhDOS is shown in 
Fig.0 

The q—0 vibration modes in a supercell include 
the "true" zone-center phonon (relative to the basic 
zincblende cell) , as well as zone-boundary phonons back- 
folded to a smaller BZ of the superstructure. It is rel- 
atively easy to distinguish them by their eigenvectors; 



certain modes are labeled in Fig. [7| Most obviously, the 
acoustical branch hits the BZ boundary at 87 cm~^ (Zu- 
rich system) or 140 cm^^ (Be-rich system). As ranges of 
frequencies of Be-related and Zn-related optical modes 
are well separated, due to a large difference in masses, it 
is straightforward to recognize the BeSe-type F-TO mode 
at 494 cm~^ (lower panel of Fig. [7|), the backfold of the 
optical branch at the zone boundary at 450 cm~^, and 
the backfold of the zone-boundary longitudinal phonon 
at 560 cm^"'^ - all in good agreement with the linear- 
response phonon dispersions calculated for BeSe by Tu- 
tuncu et al^. The zone-center LO mode is not present, 
because we do not include macroscopic electric field in 
crystal in our calculation. In the Zn-rich crystal, the 
Be-related TO mode (at 487 cm~^) behaves like an im- 
purity mode and exhibits no dispersion. Whereas the 
Zn-related mode is similarly dispersionless (at 260 cm^^) 
in the Be-rich system, in BeZn3Se4 the ZnSe-related TO 
branch shows dispersion, manifested by backfolding of 
zone-boundary phonon and hence additional structure at 
200-240 cm~^, in agree ment with experiments and cal- 
culations of Refs. I2al30i 

Even as these 8-atom supercells may be too simplistic 
to imitate short-range order effects in real quasibinary 
alloy, they exhibit important features which also persist 
in more sophisticated representative structures: a dif- 
ference between dispersionless modes due to "isolated" 
impurities and dispersing modes due to continuous con- 
nected chains; a more pronounced dispersion in BeSe- 
related modes than in much softer ZnSe-related ones. 

We turn next to the discussion of phonons in a larger 
Be6Zn26Se32 supercell. There arc now more lines in 
the spectrum, and, as the Brillouin zone (BZ) becomes 
smaller, the F modes of the supercell offer a sampling 
over more non-zone-center modes in a zincblende lattice. 
This leads in more features and more "filled" form of the 
PhDOS (Fig. El. 

The vibration spectra of both pure systems ZnSe and 
BeSe are shaped by broad acoustical branches and more 
narrow optical ones. In ZnSe they nearly overlap, re- 
sulting in a characteristic two-peaked structure of the 
PhDOS - see, e.g., Hennion et alm^. Our low- frequency 
(ZnSe-related, up to 300 cm~^) part of the supercell Ph- 
DOS provides a fair agreement with this known result 
for pure ZnSe. The spectrum of BeSe has a large separa- 
tion between its acoustical and optical parts. The former 
covers the whole region of the phonon spectrum of ZnSe. 
Therefore, in a mixed crystal, the partial PhDOS of Se 
atoms which have four Zn neighbors essentially repeats 
the shape of the Zn PhDOS, whereas those Se atoms with 
one or two Be neighbors participate in vibrations at 100 - 
200 cm^^, i.e., in the "pseudogap" region of ZnSe. These 
modes have mostly acoustical character. 

We consider now more attentively the high-frequency 
optical modes with high participation of Be vibrations. 
The corresponding PhDOS (Fig. |H1 top panel) is sepa- 
rated into contributions from (two) "isolated Be" atoms 
and those (four) in the infinite chain. The "isolated Be" 
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FIG. 8: Phonon density of states for the Be6Zn26Se32 su- 
percell, resolved over different groups of atoms, calculated for 
g=0 of the supercell and broadened by 10 cm~^. The vertical 
scaling for different groups is arbitrary. 



contribution is, expectedly, a single narrow line, coming 
about from several nearly degenerate vibration frequen- 
cies. The vibrations of Be in the chains are more diver- 
sified: the characteristic frequencies can be found below, 
within and above the "isolated Be" peak. The most re- 
markable is a clear split off at a soft side, which makes 
a distinguished peak in the vibration spectrum, similar 
to that experimentally observed in samples with a pre- 
sumed percolation on the Be sublattic o^i^^ . The previ- 
ous discussion on the distribution of bondlengths, and 
the force constants' dependence on the latter, helps us to 
understand this behavior. The bonds between the chain 
Be atoms and out-of-chain Se are shorter that those in- 
volving isolated Be, and corresponding force constants 
are larger. On the contrary, the Be-Se bonds along the 
chain are less contracted; the longer bond length implies 
smaller force constant hence lower vibration frequency. 
Our explanation of the experimentally observed anoma- 
lous (split-off) Be mode is, therefore, the following. As 
the concentration of Be in ZnSe reaches the percolation 
limit and infinite chains are built, the Be-Se bonds in 
these chains are forced to become longer that around 
isolated Be impurities, or in short chain fragments. This 
immediately softens the vibration modes which involve 



such extended bonds. 

In order to extend this interpretation over experimen- 
tally measured Raman lines, which are known to scan 
primarily the zone-center phonons, one needs yet to ex- 
tract a wavevector-resolved information about calculated 
vibrations. As the supercell size increased, the effec- 
tive sampling over the Brillouin zone (BZ) of the basic 
zincblende lattice became more fine. In particular, one 
has not only X, K and L zone-boundary phonons sam- 
pled, as was the case for the 8-atom cells, but also q points 
halfway to them from the zone center. We mentioned al- 
ready that the lack of translation symmetry in the mixed 
BeSe-ZnSe supercell does not allow to Fourier transform 
the dynamical matrix and arrive at neat phonon disper- 
sion curves. However, one can project the dispersion pat- 
terns, corresponding to different modes of vibration, onto 
the plane wave with different q values and in such way 
define q-resolved contributions to the phonon DOS: 

hico,q) = ^5]|A^(c.)exp(iqRJ|^ . (4) 

Since eigenvectors correspond to g=0 of the supercell and 
hence are real, the phase of the plane wave has no im- 
portance. One can view /i<(ci;,q) as spectral function, 
which in principle would show how the phonon dispersion 
bands of pure constituents get overlapped, distorted and 
smeared in a mixed crystal. Unfortunately, the results 
obtained for our - yet relatively small - 64-atom super- 
cell do not allow to see anything resembling a continuous 
displacement of intensity maxima with q'^^. Nonetheless, 
we present in Fig.|51the projected PhDOS for three q val- 
ues along the (001) direction. The ZnSe-type vibrations 
make clearly visible the displacement of two acoustical 
bands, known, e.g., from the measured phonon disper- 
sion curves^, as the projection wavevector changes from 
the zone-center via midpoint to the zone boundary. 

For the Be-related modes, because of small amount of 
Be atoms in the supercell, the dispersion is generally less 
pronounced, yet noticeable in the (001) direction, which 
is the general direction of the chains. Being almost non- 
interacting linear defects, the chains do not give rise to 
noticeable dispersion in perpendicular directions. The 
three acoustical modes, removed in the previous PhDOS 
figures, are retained in Fig. |51 in order to show how the 
spectral weight of the uj=0 modes disappears for q values 
away from zero. A notable difference is the PhDOS of iso- 
lated Be, which is nearly identical for q—Q and g=(001). 
This is because the supercell contains two isolated Be 
atoms. They interact only weakly but create a super- 
structure that effectively halves the BZ in the z direction. 

The g-resolving of the PhDOS singles out the modes 
which are likely to dominate the Raman spectra. The 
separation into q—0 and non-g=0 modes is however not 
very clear-cut - in part due to physical reasons (mixed 
crystal, loss of periodicity) and in part due to techni- 
cal limitations (small size of supercell and hence non- 
vanishing q=0 projection of many modes due to numer- 
ical noise). One can expect better contrast in selecting 
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FIG. 9: Phonon spectral function in the Be6Zn26Se32 super- 
cell, corresponding to three values of q of the basic zincblende 
lattice, broadened by 10 cm~^. 



q=0 modes in a larger supercell, due to better sampling. 
Nevertheless, it is clear that the separation into "isolated 
Be" and "Be in the chains" modes occurs not only on 
the average over the BZ, but is actually even more pro- 
nounced in the zone-center modes. Therefore, our expla- 
nation of anomalous Raman lines will hold. 

Most of the modes are of very mixed character, and 
involve stretching and bending of bonds in a complicated 
manner. However, we selected several ones, involving 
the -Be-Se- chains, whose vibration pattern can be de- 
scribed with relative simplicity. Two of these modes (441 
and 452 cm~^) lie in the "anomalous" Raman line, where 
the vibrations of isolated Be atoms are negligible. Three 
other selected modes, at 485, 495 and 502 cm"-'^, involve 
simultaneous vibration of all Be atoms, in the chains as 
well as isolated ones. Fig. ^] gives a snapshot of these 
five modes. Displacement patterns in some other modes 
are shown in Fig. 5 of Ref. i^. The vibrations of pre- 
dominantly "impurity Be" character, leaving the chain 
Be atoms silent, occur at nearly 490 cm~^, in between 
the third and the fourth of depicted modes. 

The most general observation is that the "anomalous" 
modes with reduced frequencies involve nearly in-plane 
displacement of Be with its own equilibrium position and 
its two Se neighbors. Such motion stretches or shortens 
the Be-Se bonds but roughly preserves the Be-Se-Be an- 



gles. The 441 cm~^ mode is a wave of correlated Be dis- 
placement along the chain, periodic with the chain step. 
This is a predominantly F mode. The 452 cm~^ mode has 
similarly "longitudinal" character but the doubled trans- 
lation length: two Be atoms approach, or depart from, 
their common Se neighbor symmetrically. Consequently 
this is not a clear zone-center vibration, however due to 
the curvature of the chain and a complexity of detailed 
displacement pattern it provides a q=0 contribution in 
the analysis according to Eq. Q). 

The harder modes are related to out-of-plane Be vibra- 
tions and hence bending of covalent Be-Se bonds. Some- 
how simplifying, in the "anomalous modes" there are pri- 
marily central Be-Se forces at work, so the frequency is 
lower than that involving an isolated Be impurity - in ac- 
cord with longer anion-cation distance along the chain, 
as illustrated by Fig. El In the latter case, even as in- 
teratomic distances remain the same, the off-center con- 
tributions due to the bond bending effectively increase 
the force constants and harden the resulting frequency 
beyond those for an isolated Be. 



VII. CONCLUSIONS 

We simulated the equilibrium structure and lattice dy- 
namics from first principles in a series of supercells simu- 
lating mixed (Be,Zn)Se crystals. We found a nearly linear 
dependence of equilibrium lattice constant with compo- 
sition, but a relative independence of (average) Be-Se 
and Zn-Se distances on concentration. When incorpo- 
rated in a mixed crystal, the individual interatomic dis- 
tances develop a certain scattering around their mean 
values, maintaining however a clear separation between 
the Be-Se and Zn-Se bond lengths. Specifically, an iso- 
lated Be atom may much more efficiently shorten the 
bonds to an Se neighbor whose other three bonds are 
Zn-terminated, than to an Se atom shared by other Be 
neighbors. In particular, a continuous -Be-Se- chain is 
caracterized by longer link lengths than the average Be- 
Se distance, for a given concentration. We demonstrated 
further that this has effect on force constants and vibra- 
tion frequencies. As follows from our direct calculation 
of interatomic force constants, the principal values of the 
latter, taken as diagonal elements of a 3x3 matrix re- 
lating two selected atoms, show a linear decrease as a 
function of cation-anion distance between nearest neigh- 
bors. The calculations manifest a primarily ionic char- 
acter of the Zn-Se interaction and a remarkable level of 
covalency in the Be-Se bonds. A decrease of force con- 
stants with larger bond distance is the microscopic origin 
of the formation of a softened "anomalous" mode, which 
appears on the low-frequency side of BeSe-related TO 
mode. Our wavevector analyzis of vibration patterns in 
the supercell indicates that the "anomalous" mode has a 
strong zone-center contribution, that is consistent with 
its clear experimental observation in the Raman spectra. 
We confirmed therefore on the microscopic level the early 
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FIG. 10: Vibration patterns of five selected modes with substantial contribution of Be atoms; the relative displacements of 
the latter are shown by arrows. Only some atoms of the Be6Zn26Se32 supercell (Fig.|5] right panel) are shown. The figure was 
created with the XCrySDen softwareS™. The first two modes affect exclusively the continuous -Be-Se- chains; the last three 
contain also comparable contribution of impurity Be atoms. See text for details. 

presumed relation between the onset of Be percolation de Physique. Useful comments by N. E. Christensen, 
threshold on the cation sublattice and the appearance of discussions with V. Vikhnin and travel support from the 
"anomalous" mode, including specific vibrations of con- NATO project CLG 980378 are greatly appreciated, 
tinuous -Be-Se- chains transversing the crystal. 
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